/*
 * Voronoi.cpp
 *
 *  Created on: Jun 20, 2011
 *      Author: marchi
 */

#include "Voronoi.h"
int Voronoi::nresid=0;
int Voronoi::nr=0;
int Voronoi::nc=0;
float Voronoi::time=0.0;
string * Voronoi::label=NULL;
int * Voronoi::types=NULL;
container_periodic_poly * Voronoi::MyconP=NULL;
container_poly * Voronoi::MyconNP=NULL;
wall_initial_shape * Voronoi::wis0=NULL;
bool Voronoi::Spheres=false;
vector<vector<int> > dummy;


struct op_SumMinus{
	struct duo{
		int n;
		double a;
	};
	struct isNotMinus{bool operator()(const duo & x){return x.n>0;};};
	void operator()(vector<int> & In,vector<double> & A){
		vector<duo> tmp=vector<duo>(In.size());
		for(size_t o=0;o < tmp.size();o++){
			tmp[o].n=In[o];
			tmp[o].a=A[o];
		}
		vector<duo>::iterator divide=std::stable_partition(tmp.begin(),tmp.end(),isNotMinus());
		vector<duo>::iterator it=divide;
		size_t Newdim=std::distance(tmp.begin(),divide)+1;
		double TArea=0.0;
		for(;it != tmp.end(); it++) {
			TArea+=it->a;
		}
		In.resize(Newdim);
		A.resize(Newdim);
		it=tmp.begin();
		for(size_t o=0;o<Newdim-1;o++){
			In[o]=tmp[o].n;
			A[o]=tmp[o].a;
		}
		In[Newdim]=tmp[Newdim-1].n;
		A[Newdim-1]=TArea;
	}
};
struct tArea{
	int n;
	double a;
	tArea(): a(0), n(0){};
	tArea(int m){n=m;a=0.0;}
	tArea(int m, double b){n=m;a=b;}
	virtual ~tArea(){};
};

struct tAcomp{
	bool operator()(const tArea & x, const tArea & y) const {return x.n < y.n;};
};
struct isMinus{
	bool operator()(const tArea & x) const {return x.n < 0;};
};

int wall_initial_shape::count=0;
using namespace Topol_NS;
Voronoi::Voronoi(Topol & y, bool bH): TypesName(NULL),Vol(NULL),porder(NULL),
		Neighs(NULL),Surface(NULL){
	nr=y.Size();
	Vol=new double [nr];
	Centroids=vector<Dvect>(nr);
	Neighs=new vector<int> [nr];
	Surface=new vector<double> [nr];
	label=new string [nr];
	types=new int [nr];
	CIndex=y.gCIndex();
	for(size_t o=0;o< nr;o++){
		Vol[o]=0.0;
		label[o]=y.AtomInfo(o);
		if(bH || label[o].compare(0,1,h,0,1)) {
			cindex.push_back(o);
		}
	}
	Atominfo=vector<int>(nr);
	for(size_t o=0;o<CIndex.size();o++)
		for(size_t p=0;p<CIndex[o].size();p++){
			int i=CIndex[o][p];
			Atominfo[i]=o;
		}

	restype_s=y.getrestype_s();
	restype_i=y.getrestype_i();
	ResType=y.getResType();
	nc=y.getmaxt();
	nresid=y.ResSize();
	area.Allocate(nresid,nc+1);
	nearest.Allocate(nresid,nc+1);
	Vols.Allocate(nresid);
	Cxcm=vector<Dvect>(nresid);
	TypesName=new string [nc];
	const string * tname=y.getTypeNames();
	for(int o=0;o<nc;o++){
		TypesName[o]=tname[o];
		const vector<int> & cidx=y.getCidx(o);
		for(size_t ia=0;ia<cidx.size();ia++){
			types[cidx[ia]]=o;
		}
	}
	for(int o=0;o<nresid;o++)
		RealResidue.push_back(o);
	resinfo=vector<string>(nresid);
	for(int o=0;o<nresid;o++)
		resinfo[o]=y.ResInfo(o);
}

Voronoi::Voronoi(): TypesName(NULL),Vol(NULL),porder(NULL),
		Neighs(NULL),Surface(NULL) {
}
Voronoi::Voronoi(int n, const vector<string> & s, bool bH):  TypesName(NULL),
		Vol(NULL),porder(NULL),Neighs(NULL),Surface(NULL){
	nr=n;
	Vol=new double [nr];
	Centroids=vector<Dvect>(nr);
	Neighs=new vector<int> [nr];
	Surface=new vector<double> [nr];
	label=new string [nr];
	types=new int [nr];
	for(int i=0;i<nr;i++){
		Vol[i]=0.0;
		label[i]=s[i];
		if(bH || label[i].compare(0,1,h,0,1)) cindex.push_back(i);
	}
}
Voronoi::Voronoi(int n, char *** s, bool bH): TypesName(NULL),Vol(NULL),porder(NULL),
		Neighs(NULL),Surface(NULL), CIndex(dummy){
	nr=n;
	Vol=new double [nr];
	Centroids=vector<Dvect>(nr);
	Neighs=new vector<int> [nr];
	Surface=new vector<double> [nr];
	label=new string [nr];
	types=new int [nr];
	for(int i=0;i<nr;i++){
		Vol[i]=0.0;
		label[i]=*s[i];
		if(bH || label[i].compare(0,1,h,0,1)) cindex.push_back(i);
	}
}
inline bool compare(Atoms::plane x,Atoms::plane y){
	return x.xc[3] < y.xc[3];
}

void Voronoi::setContacts(TopolPDB & y){
	TrueRes=vector<string>(y.Size());
	try{
		stringstream ss1,ss2;
		ss1<<nr;
		ss2<<y.Size();
		if(nr != y.Size()) throw string(" Number of atoms is "+ ss1.str() +
				" while from reading PDB " + ss2.str() + " atoms were found. Abort.");
	} catch(const string & s){
		cout << s << endl;
		exit(1);
	}
	for(size_t o=0; o< y.Size(); o++){
		string sub1=y[o].first.substr(11,5);
		string sub2=y[o].first.substr(17,3);
		string sub3=y[o].first.substr(21,1);
		string sub4=y[o].first.substr(23,4);
		sub1.erase(remove_if(sub1.begin(),sub1.end(),::isspace),sub1.end());
		sub2.erase(remove_if(sub2.begin(),sub2.end(),::isspace),sub2.end());
		sub3.erase(remove_if(sub3.begin(),sub3.end(),::isspace),sub3.end());
		sub4.erase(remove_if(sub4.begin(),sub4.end(),::isspace),sub4.end());
		stringstream ss;
		ss<< setw(3) << left << fixed << sub1 << " ";
		ss<< setw(3) << right << fixed << sub2 << "-";
		ss<< setw(5) << left << fixed << sub3 + sub4;;
		TrueRes[o]=ss.str();
	}
}
void Voronoi::StartP(float frame, Atoms & atm){
	time=frame;
	const int nx=NNN,ny=NNN,nz=NNN;
	for(int o=0;o<DIM;o++)
		for(int p=0;p<DIM;p++){
			co[o][p]=atm.getMt().getCO()[o][p];
			oc[o][p]=atm.getMt().getOC()[o][p];
		}
	FAtoms * atms=dynamic_cast<FAtoms *> (&atm);

	double bx=co[0][0], bxy=co[0][1], by=co[1][1],bxz=co[0][2],byz=co[1][2],bz=co[2][2];

	vector<double> vertx;
	if(MyconP) delete MyconP;
	if(porder) delete porder;
	MyconP=new container_periodic_poly(bx,bxy,by,bxz,byz,bz,nx,ny,nz,8);
	porder=new particle_order;
	for(unsigned int o=0;o<cindex.size();o++){
		double x=atm[cindex[o]][XX];
		double y=atm[cindex[o]][YY];
		double z=atm[cindex[o]][ZZ];
		double r=atms->getrd(cindex[o]);
		MyconP->put(*porder,o,x,y,z,r);
	}
}
void Voronoi::StartNP(float frame, Atoms & atm){
	time=frame;
	if(!wis0) wis0=new wall_initial_shape ;

	const int nx=NNN,ny=NNN,nz=NNN;
	double sum=0.0;
	for(int o=0;o<DIM;o++)
		for(int p=0;p<DIM;p++){
			co[o][p]=atm.getMt().getCO()[o][p];
			oc[o][p]=atm.getMt().getOC()[o][p];
			if(o<p) sum+=abs(co[o][p]);
		}
	try{
		if(abs(sum) >SMALL) throw string("Can only work with orthogonal axis. Abort! ");
	} catch(const string & s){
		cout << s <<endl;
		exit(1);
	}
	FAtoms * atms=dynamic_cast<FAtoms *> (&atm);

	double bx=0.5*co[XX][XX], by=0.5*co[YY][YY],bz=0.5*co[ZZ][ZZ];

	vector<double> vertx;
	if(MyconNP) delete MyconNP;
	if(porder) delete porder;
	MyconNP=new container_poly(-bx,bx,-by,by,-bz,bz,nx,ny,nz,true,true,true,8);
	porder=new particle_order;

	MyconNP->add_wall(*wis0);
	for(unsigned int o=0;o<cindex.size();o++){
		double x=atm[cindex[o]][XX];
		double y=atm[cindex[o]][YY];
		double z=atm[cindex[o]][ZZ];
		double r=atms->getrd(cindex[o]);
		MyconNP->put(*porder,o,x,y,z,r);
	}
}
void Voronoi::StartNP(Atoms & atm){
	time=0;
	if(!wis0) wis0=new wall_initial_shape ;
	const int nx=NNN,ny=NNN,nz=NNN;
	for(int o=0;o<DIM;o++)
		for(int p=0;p<DIM;p++){
			co[o][p]=atm.getMt().getCO()[o][p];
			oc[o][p]=atm.getMt().getOC()[o][p];
		}
	FAtoms * atms=dynamic_cast<FAtoms *> (&atm);

	double bx=100.0, by=100.0,bz=100.0;

	vector<double> vertx;
	if(MyconNP) delete MyconNP;
	if(porder) delete porder;
	MyconNP=new container_poly(-bx,bx,-by,by,-bz,bz,nx,ny,nz,false,false,false,8);
	porder=new particle_order;

	MyconNP->add_wall(*wis0);
	for(unsigned int o=0;o<cindex.size();o++){

		double x=atm[cindex[o]][XX];
		double y=atm[cindex[o]][YY];
		double z=atm[cindex[o]][ZZ];
		double r=atms->getrd(cindex[o]);
		MyconNP->put(*porder,o,x,y,z,r);
	}
}
void Voronoi::getDelaunay(Cavities::isTypes & y){}

void Voronoi::testVol(){
	double VolCell=Metric(co).getVol();
	double VorVol=0.0;
	for(unsigned int n=0;n<cindex.size();n++)
		VorVol+=Vol[cindex[n]];
	cout << setw(10) << setprecision(2) << scientific
			<< "Volume error is " << 1000.0*(VorVol-VolCell)
			<< " A^3 over " << setprecision(4) << fixed << 1000.0*VolCell << " A^3 "<< endl;

}
void Voronoi::gather(vector<int> & it){
	vector<int> tmp=it;
	for(unsigned int o=0;o<it.size();o++) {
		if(it[o]> -1) tmp[o]=cindex[it[o]];
	}
	it=tmp;
	tmp.clear();
}
void Voronoi::setResTypes(Topol_NS::Topol & MyTop){
	restype_s=MyTop.getrestype_s();
	restype_i=MyTop.getrestype_i();
	ResType=MyTop.getResType();
}

void Voronoi::setTypes(int nc0, const AtomIndex cidx[], const string * tname){
	nc=nc0;
	nresid=ResidueCM::Size();
	area.Allocate(nresid,nc+1);
	nearest.Allocate(nresid,nc+1);
	Vols.Allocate(nresid);
	Cxcm=vector<Dvect>(nresid);
	TypesName=new string [nc];
	for(int o=0;o<nc;o++){
		TypesName[o]=tname[o];
		for(int ia=0;ia<cidx[o].getN();ia++){
			types[cidx[o][ia]]=o;
		}
	}
	for(int o=0;o<ResidueCM::getiDSize();o++)
		RealResidue.push_back(ResidueCM::getiD(o));
}
void Voronoi::testit(){
	double vol0;
	c_loop_order vl(*MyconNP,*porder);
	voronoicell_neighbor c;
	int ia=0;
	if(vl.start())
		do {
			if(MyconNP->compute_cell(c,vl)) {
				}
				ia++;
				cout << c.volume() << endl;
		} while(vl.inc());

	exit(1);

}

void Voronoi::getData(vector<Topol_NS::MySubunit> & MySub){
	if(Phobic.size() == 0){
		vector<int> tmp(MySub.size());
		vector<string> tmp2(MySub.size());
	 	for(size_t i=0;i<MySub.size();i++){
	 		tmp[i]=MySub[i].getN();
	 		tmp2[i]=MySub[i].getC();
	 	}
	 	vector<string>::iterator it2;
	 	it2=unique(tmp2.begin(),tmp2.end());
	 	tmp2.resize(distance(tmp2.begin(),it2));
	 	for(size_t i=0;i<tmp2.size();i++)
	 		unique_domains.push_back(tmp2[i]);


	 	size_t m=*max_element(tmp.begin(),tmp.end())+1;
	 	Phobic=vector<vector<double> >(m);
	 	Philic=vector<vector<double> >(m);
	 	Mixt=vector<vector<double> >(m);
	 	for(size_t n=0;n<m;n++){
	 		Phobic[n]=vector<double>(m);
	 		Philic[n]=vector<double>(m);
	 		Mixt[n]=vector<double>(m);
	 	}
	}
	size_t m=Phobic.size();
	for(size_t n=0;n<m;n++){
		for(size_t p=0;p<m;p++){
			Phobic[n][p]=0.0;
			Philic[n][p]=0.0;
			Mixt[n][p]=0.0;
	 	}
	}


	for(unsigned int n=0;n<cindex.size();n++){
		Vol[cindex[n]]=0.0;
		Neighs[cindex[n]].clear();
		Surface[cindex[n]].clear();
	}


	double AreaA=0.0;
	vector<int> nei;
	vector<int> order0;
	vector<double> area0;
	voronoicell_neighbor c0,c;
	int ia=0;
	double vol0;
	area=0.0;
	nearest=0;
	if(MyconP) {
		c_loop_order_periodic vl(*MyconP,*porder);
		if(vl.start())
			do {
				if(MyconP->compute_cell(c,vl)) {
					nei.clear();
					area0.clear();
					vol0=c.volume();
					c.neighbors(nei);

					c.face_areas(area0);
					double cx,cy,cz;
					c.centroid(cx,cy,cz);
					double x,y,z;
					int yu;
					double rt;
					vl.pos(yu,x,y,z,rt);
					gather(nei);
					Centroids[cindex[ia]]=Dvect(x+cx,y+cy,z+cz);
					Vol[cindex[ia]]=vol0;
					Neighs[cindex[ia]]=nei;
					Surface[cindex[ia]]=area0;
					Topol_NS::MySubunit a1=MySub[cindex[ia]];
					for(unsigned int p=0;p<nei.size();p++){
						Topol_NS::MySubunit a2=MySub[nei[p]];
						int n1=a1.getN();
						int n2=a2.getN();
						if(n1 < n2) {
							if(a2.isHydrophobic(a1)) {
								Phobic[n1][n2]+=area0[p];
							} else if(a2.isHydrophilic(a1)) {
								Philic[n1][n2]+=area0[p];
							} else if(a2.isMixed(a1)) {
								Mixt[n1][n2]+=area0[p];
							}
						}
					}
				}
			ia++;
		} while(vl.inc());
	} else if(MyconNP) {
		double volumetot=0.0;
		c_loop_order vl(*MyconNP,*porder);
		if(vl.start())
			do {
				int yu;
				double x,y,z,rt;
				vl.pos(yu,x,y,z,rt);
				wis0->setR0(rt);
				//wis0->setR0(0.45);
				if(MyconNP->compute_cell(c,vl)) {
					nei.clear();
					area0.clear();
					vol0=c.volume();
					c.neighbors(nei);
					c.face_areas(area0);
					op_SumMinus()(nei,area0);
					double cx,cy,cz;
					c.centroid(cx,cy,cz);
					gather(nei);
					Centroids[cindex[ia]]=Dvect(x+cx,y+cy,z+cz);
					Vol[cindex[ia]]=vol0;
					Neighs[cindex[ia]]=nei;
					Surface[cindex[ia]]=area0;
					Topol_NS::MySubunit a1=MySub[cindex[ia]];
					for(unsigned int p=0;p<nei.size();p++){
						if(nei[p] <0) continue;
						Topol_NS::MySubunit a2=MySub[nei[p]];
						int n1=a1.getN();
						int n2=a2.getN();
						if(n1 < n2) {
							if(a2.isHydrophobic(a1)) {
								Phobic[n1][n2]+=area0[p];
							} else if(a2.isHydrophilic(a1)) {
								Philic[n1][n2]+=area0[p];
							} else if(a2.isMixed(a1)) {
								Mixt[n1][n2]+=area0[p];
							}
						}
					}
				}
			ia++;
			volumetot+=vol0;
		} while(vl.inc());
	}

	nei.clear();
	area0.clear();
	vector<int> dummy;
	for(int o=0;o<nc;o++)
		for(int p=0;p<nc;p++)
			mapIntRes[TypesName[o]][TypesName[p]]=0.0;
	TrueNei=vector<vector<int> >(nr);
	for(int o=0;o<nresid;o++) {
		vector<int> & cindex_o=dummy;
		if(!ResidueCM::isNULL()) cindex_o=ResidueCM::getind(o);
		if(!CIndex.empty()) cindex_o=CIndex[o];
		double sum_v=0.0;
		Dvect sum_cm=Dvect(0.0,0.0,0.0);
		vector<tArea> cdx(cindex_o.begin(),cindex_o.end());
		sort(cdx.begin(),cdx.end(),tAcomp());
		for(unsigned int ia=0;ia<cdx.size();ia++){
			int i=cdx[ia].n;
			if(Vol[i] != 0.0) {
				int na=cdx[ia].n;
				for(unsigned int ib=ia+1;ib<cdx.size();ib++){
					int j=cdx[ib].n;
					if(Vol[j] != 0.0) {
						int nb=cdx[ib].n;
						Centroids[nb].Displ(Centroids[na],co,oc);
					}
				}
			}
		}
		int o_type=getTypesRes(o);
		vector<tArea> nei0,nei;

		for(unsigned int ia=0;ia<cdx.size();ia++){
			nei0.clear();nei.clear();
			int i=cdx[ia].n;
			sum_v+=Vol[i];
			sum_cm+=Centroids[i]*Vol[i];
			for(unsigned int p=0;p<Neighs[i].size();p++){
				nei0.push_back(tArea(Neighs[i][p],Surface[i][p]));
			}
			sort(nei0.begin(),nei0.end(),tAcomp());
			set_difference(nei0.begin(),nei0.end(),cdx.begin(),cdx.end(),back_inserter(nei),tAcomp());
			vector<tArea>::iterator it=nei.begin();
			for(;it != nei.end(); ++it){
				//if(getTypes(i) == getTypes(it->n) && i <= it->n) continue;
				if(it->n > -1) {
					TrueNei[i].push_back(it->n);
					area[o][getTypes(it->n)]+=it->a;
					nearest[o][getTypes(it->n)]++;
				} else {
					area[o][nc]+=it->a;
					nearest[o][nc]++;
				}
			}

		}
		string rinfo=resinfo[o];
		for(int p=0;p<nc;p++){
			mapIntRes[rinfo][TypesName[p]]+=area[o][p];
		}
		Vols[o]=sum_v;
		sum_cm=sum_cm/sum_v;
		Cxcm[o]=sum_cm;
	}
	/*
	for(size_t o=0;o< nresid;o++)
		if(area[o][nc]> 5.0e-3) cout << o+1 << " " ;
	cout << endl;
	cout << " core " <<endl;
	for(size_t o=0;o< nresid;o++)
		if(area[o][nc] < 5.0e-3) cout << o+1 << " " << area[o][nc]<<endl;

	exit(1);
	*/
	for(int o=0;o<nc;o++)
		mapIntRes[TypesName[o]][TypesName[o]]*=0.5;
}
void Voronoi::WriteIt(std::ofstream & fout){
	fout << "######>> At step No. " << setw(10) << setprecision(2) << fixed<< time << endl;
	int po=0;
	for(int o=0;o<nresid;o++) {
		if(!VoronoiPrint::bPrintVols)continue;
		if(restype_i[getTypesRes(o)] != VoronoiPrint::pGroup && VoronoiPrint::pGroup != -1) continue;
		double a=Vols[o]*1000.0;
		string l=(ResidueCM::isNULL())?resinfo[o]:ResidueCM::getlabels(o);
		int rres=RealResidue[o]+1;
		(po%5)?fout << setw(10) << setprecision(4) << fixed<<  a << ' ' << setw(5) << fixed << l << ' ' << setw(5) << fixed << rres << ' ':
		  fout << endl << "%$VolRes " << setw(10) << setprecision(4) << fixed << a << ' ' << setw(5) << fixed << l  <<' ' << setw(5) << fixed << rres << ' ';
		po++;
	}
	fout << endl;
	array2<double> interface;
	array2<int> innei;
	interface.Allocate(nc,nc+1);
	innei.Allocate(nc,nc+1);
	interface=0.0;
	innei=0;
	array1<double> VolSel;
	VolSel.Allocate(nc);
	VolSel=0.0;
	for(int o=0;o<nresid ;o++){
		int o_type=getTypesRes(o);
		VolSel[o_type]+=Vols[o]*1000.0;
		for(int p=0;p<nc+1;p++) {
			double a=area[o][p]*100.0;
			interface[o_type][p]+=a;
			innei[o_type][p]+=nearest[o][p];
		}

	}

	if(VoronoiPrint::bPrintDomain){
		size_t m=Phobic.size();
		fout << "# Hydrophobic Surface Area:\n";
		for(size_t n=0;n<m;n++){
			Phobic[n][n]=0.0;
			Philic[n][n]=0.0;
			Mixt[n][n]=0.0;
			for(size_t p=n+1;p<m;p++){
				Phobic[p][n]=Phobic[n][p];
				Philic[p][n]=Philic[n][p];
				Mixt[p][n]=Mixt[n][p];
			}
		}
		for(size_t n=0;n<m;n++){
			string s1=unique_domains[n];
			for(size_t p=0;p<m;p++){
				string s2=unique_domains[p];
				if(p == 0) fout << " Phobic ";
				fout << fixed <<" " << s1 << " -- " << s2 << setprecision(3) << fixed << " = " << setw(10) << Phobic[n][p]*100.0 << "   ";
		 	}
			fout << "\n";
		}
		fout << "\n";
		fout << "# Hydrophilic Surface Area:\n";
		for(size_t n=0;n<m;n++){
			string s1=unique_domains[n];
			for(size_t p=0;p<m;p++){
				string s2=unique_domains[p];
				if(p == 0) fout << " Philic ";
				fout << fixed <<" " << s1 << " -- " << s2 << setprecision(3) << fixed << " = " << setw(10) << Philic[n][p]*100.0 << "   ";
		 	}
			fout << "\n";
		}
		fout << "\n";
		fout << "# Mixed Surface Area:\n";
		for(size_t n=0;n<m;n++){
			string s1=unique_domains[n];
			for(size_t p=0;p<m;p++){
				string s2=unique_domains[p];
				if(p == 0) fout << " Mixed ";
				fout << fixed <<" " << s1 << " -- " << s2 << setprecision(3) << fixed << " = " << setw(10) << Mixt[n][p]*100.0 << "   ";
		 	}
			fout << "\n";
		}
		fout << "\n Total Areas:  " << endl;
		vector<double> Tphobic(m);
		vector<double> Tphilic(m);
		vector<double> Tmixt(m);
		for(size_t n=0;n<m;n++){
			Tphobic[n]=0.0;
			Tphilic[n]=0.0;
			Tmixt[n]=0.0;
			string s1=unique_domains[n];
			for(size_t p=0;p<m;p++){
				Tphobic[n]+=Phobic[n][p];
				Tphilic[n]+=Philic[n][p];
				Tmixt[n]+=Mixt[n][p];
			}
			double Tot=Tphobic[n]+Tphilic[n]+Tmixt[n];
			fout << fixed <<" " << s1 << "     :  ";
			fout << " Phobic " << setprecision(3) << fixed << " -> " << setw(10) << Tphobic[n]*100.0 << "   ";
			fout << " Philic " << setprecision(3) << fixed << " -> " << setw(10) << Tphilic[n]*100.0 << "   ";
			fout << " Mixed "  << setprecision(3) << fixed << " -> " << setw(10) << Tmixt[n]*100.0 << "   ";
			fout << " Total "  << setprecision(3) << fixed << " = " << setw(10) << Tot*100.0 << " \n";
		}
	}

	if(VoronoiPrint::bPrintAreas ){
		vector<double> area_r=vector<double>(Topol_NS::ResidueTypes::getTotTypes(),0.0);
		vector<int> nei_r=vector<int>(Topol_NS::ResidueTypes::getTotTypes(),0);
		fout << "# Area format:        ";
		for(int o=0;o<ResType.size();o++)
			fout << setw(1) << fixed << " areatype[" << Topol_NS::ResidueTypes::getType(ResType[o]) << "(" << ResType[o] << ")"<< "] ";
		fout << endl;
		fout  << "%$AreaRes " ;
		int po=0;
		double AreaNotCovered=0.0;
		for(int o=0;o<nresid ;o++){
			string l=(ResidueCM::isNULL())?resinfo[o]:ResidueCM::getlabels(o);
			int o_type=getTypesRes(o);
			if(VoronoiPrint::pGroup != -1 && restype_i[o_type] != VoronoiPrint::pGroup) continue;
			fout  << right << setw(5)<< l << " " << setw(5) << fixed << RealResidue[o]+1 << ' ' << setw(3) << fixed << restype_i[o_type] << ' ' ;
			for(int p=0;p<area_r.size();p++)
				area_r[p]=0.0;
			for(int p=0;p<nc;p++) {
				int n=restype_i[p];
				area_r[restype_i[p]]+=100.0*area[o][p];

			}
			AreaNotCovered+=100.0*area[o][nc];

			for(int p=0;p<ResType.size();p++) {
				double a=area_r[ResType[p]];
				fout << setw(10) << setprecision(5) << fixed << a << ' ';
			}
			if(o+1==nresid) fout << endl;
			else if(!((po+1)%2)) fout << endl << "%$AreaRes " ;
			po++;
		}
		if(AreaNotCovered) fout << "\n%$exposed Area = " << AreaNotCovered << "\n" <<endl;
		po=0;
		fout <<endl;
		fout << "\n# Inter residue share surface:    " <<endl;;

		fout << endl << "%$IRes  ";
		for(int o=0;o<nc;o++){
			string Lab1=TypesName[o];
			for(int p=o;p<nc;p++){
				string Lab2=TypesName[p];
				double a=mapIntRes[Lab1][Lab2]*100.0;
				fout << setw(4) << fixed << left << Lab1 << "- " ;
				fout << setw(4) << fixed << left << Lab2 << " : ";
				fout << setw(9) << setprecision(3) << right << fixed << a << "  ";
				if((po+1)%3 == 0) {fout << endl << "%$IRes  ";}
				po++;
			}
		}
		fout << endl;
		fout << endl;
		for(int o=0;o<ResType.size();o++)
			fout << setw(1) << fixed << " NeigType[" << Topol_NS::ResidueTypes::getType(ResType[o]) << "(" << ResType[o] << ")"<< "] ";
		fout << endl;
		fout  << "%$NeighsRes " ;
		po=0;
		if(MyconNP){
			fout << "\n# Neighbors format:    ";
			for(int o=0;o<ResType.size();o++)
				fout << setw(1) << fixed << " NeigType[" << Topol_NS::ResidueTypes::getType(ResType[o]) << "(" << ResType[o] << ")"<< "] ";
			fout << endl;
			fout  << "%$NeighsRes " ;
			po=0;

			for(int o=0;o<nresid ;o++){
				string l=(ResidueCM::isNULL())?resinfo[o]:ResidueCM::getlabels(o);
				int o_type=getTypesRes(o);
				if(VoronoiPrint::pGroup != -1 && restype_i[o_type] != VoronoiPrint::pGroup) continue;
				fout  << right << setw(5)<< l << " " << setw(5) << fixed << RealResidue[o]+1 << ' ' << setw(3) << fixed << restype_i[o_type] << ' ' ;
				for(int p=0;p<nei_r.size();p++)
					nei_r[p]=0.0;
				for(int p=0;p<nc;p++) {
					nei_r[restype_i[p]]+=nearest[o][p];
				}
				for(int p=0;p<ResType.size();p++) {
					int a=nei_r[ResType[p]];
					fout << setw(10) << fixed << a << ' ';
				}
				if(o+1==nresid) fout << endl;
				else if(!((po+1)%2)) fout << endl << "%$NeighsRes " ;
				po++;
			}
		}
	}

	if(!VoronoiPrint::bPrintDomain) {
		int ntype=Topol_NS::ResidueTypes::getTotTypes();

		vector<double> vol_r=vector<double>(ntype,0.0);
		vector<vector<double> > inter_r=vector<vector<double> >(ntype,vol_r);
		vector<vector<int> > innei_r=vector<vector<int> >(ntype,vector<int>(ntype,0));

		fout << "\n"<<endl;
		fout << "#  Residue Types: ";
		for(int o=0;o<ResType.size();o++)
			fout << setw(1) << fixed << " areatype[" << Topol_NS::ResidueTypes::getType(ResType[o]) << "(" << ResType[o] << ")"<< "] ";
		fout << endl;

		fout << "# Volume of selection " << endl <<"#           ";
		for(int o=0;o<ResType.size();o++)
			fout << setw(1) <<fixed << " Volumetype[" << Topol_NS::ResidueTypes::getType(ResType[o]) << "(" << ResType[o] << ")"<< "] ";
		fout << endl;
		fout << "%$TotVol ";

		for(int o=0;o<nc;o++){
			int oo=restype_i[o];
			vol_r[oo]+=VolSel[o];
			for(int p=0;p<nc;p++){
				int pp=restype_i[p];
				inter_r[oo][pp]+=interface[o][p];
			}
		}
		for(int p=0;p<nc;p++){
			int pp=restype_i[p];
			inter_r[pp][pp]*=0.5;
		}
		for(int o=0;o<ResType.size();o++)
			fout << setw(21) << setprecision(4) << fixed << vol_r[ResType[o]];

		fout << "\n" <<endl;

		fout << "# Interface area of selection " << endl << "#                ";
		for(int o=0;o<ResType.size();o++)
			fout << setw(1) << fixed << " areatype[" << Topol_NS::ResidueTypes::getType(ResType[o]) << "(" << ResType[o] << ")"<< "]";
		fout << endl;
		for(int o=0;o<ResType.size();o++){
			int oo=ResType[o];
			fout << "%$AreaTot " << setw(3) << oo << "     ";
			for(int p=0;p<o;p++){
				fout <<setw(17) << setprecision(3) << fixed << " ";
			}
			for(int p=o;p<ResType.size();p++){
				int pp=ResType[p];
				fout << setw(17) << inter_r[oo][pp] ;
			}
			fout << endl;
		}
		for(int o=0;o<nc;o++){
			int oo=restype_i[o];
			vol_r[oo]+=VolSel[o];
			for(int p=0;p<nc;p++){
				int pp=restype_i[p];
				innei_r[oo][pp]+=innei[o][p];
			}
		}
		for(int p=0;p<nc;p++){
			int pp=restype_i[p];
			innei_r[pp][pp]*=0.5;
		}

		fout << "\n" <<endl;

		fout << "# Neighbors of selection=      " << endl << "#                ";
		for(int o=0;o<ResType.size();o++)
			fout << setw(1) << fixed << " NeighType[" << Topol_NS::ResidueTypes::getType(ResType[o]) << "(" << ResType[o] << ")"<< "]";
		fout << endl;
		for(int o=0;o<ResType.size();o++){
			int oo=ResType[o];
			fout << "%$TotNeigh " << setw(3) << oo << "     ";
			for(int p=0;p<o;p++){
				fout <<setw(17) << setprecision(3) << fixed << " ";
			}
			for(int p=o;p<ResType.size();p++){
				int pp=ResType[p];
				fout << setw(17) << innei_r[oo][pp] ;
			}
			fout << endl;
		}

		fout << "\n\n" <<endl;
		fout << string(90,'-') + "\n";
		fout << "# Closest contacts :      " << endl << "#                \n";
		if(!TrueRes.empty()){
			for(size_t o=0; o < cindex.size(); o++){
				int i=cindex[o];
				string sub_i=TrueRes[i].substr(0,3);
				if(restype_i[getTypes(i)] == 0){
					if(!Neighs[i].size()) continue;
					if(sub_i.find("C  ") != string::npos || sub_i.find("CA ") != string::npos || sub_i.find("N  ") != string::npos || sub_i.find("O  ") != string::npos) continue;
					int nat=0;
					for(unsigned int p=0;p<Neighs[i].size();p++){
						int pp=Neighs[i][p];
						if(pp == -1) continue;
						if(restype_i[getTypes(pp)] == 0){
							string sub_j=TrueRes[pp].substr(0,3);
							if(sub_j.find("C  ") != string::npos || sub_j.find("CA ") != string::npos || sub_j.find("N  ") != string::npos || sub_j.find("O  ") != string::npos) continue;
							if(abs(Atominfo[i]-Atominfo[pp]) > 3) nat++;
						}
					}
					if(!nat) continue;
					fout << "%$C " <<fixed << left << setw((13)) <<TrueRes[i] << " : " ;
					nat=0;
					for(unsigned int p=0;p<Neighs[i].size();p++){
						int pp=Neighs[i][p];
						if(pp == -1) continue;
						if(restype_i[getTypes(pp)] == 0){
							string sub_j=TrueRes[pp].substr(0,3);
							if(sub_j.find("C  ") != string::npos || sub_j.find("CA ") != string::npos || sub_j.find("N  ") != string::npos || sub_j.find("O  ") != string::npos) continue;
							if(abs(Atominfo[i]-Atominfo[pp]) > 3){
								if(nat && !((nat*13)%7)) fout<<"\n" + string(17,' ')+": ";
								fout << fixed << left << setw((13)) << TrueRes[pp];
								nat++;
							}
						}
					}
					fout <<endl;
				} else{
					int nat=0;
					fout << "%$C " <<fixed << left << setw((13)) <<TrueRes[i] << " : " ;
					for(unsigned int p=0;p<Neighs[i].size();p++){
						int pp=Neighs[i][p];
						if(pp == -1) continue;
						if(nat && !((nat*13)%7)) fout<<"\n" + string(17,' ')+": ";
						fout << fixed << left << setw((13)) << TrueRes[pp];
						nat++;
					}
					fout <<endl;
				}
			}
		}
	}
}


